Water Resources in Copenhagen during 20th century - 01

This script visualizes the spatial component of the data accompanying the Spring 2021 course on the City: Between Culture and Nature, taught by Mikkel Thelle and Mikkel Høghøj. The course surveys the gradual appearance of private and public bathing facilities, toilets and communal hygienic resources in the city of Copenhagen during the 20th century. By editing elements in this script, you can explore different aspects of past and present hygienic amenities in the capital of Denmark.

Before we start: data wrangling

First load the packages necessary for spatial data visualisation and analysis.

library(sf)
library(tidyverse)
library(googlesheets4)
library(leaflet)

Spatial data

Next, load your spatial data - polygons representing the suburbs of Copenhagen.

suburbs <- st_read("data/bydel.shp", options = "ENCODING=WINDOWS-1252")  # thank you, Malte, for the options argument which fixed the Danish chars

options:        ENCODING=WINDOWS-1252 
Reading layer `bydel' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\BetweenCultureAndNature\data\bydel.shp' using driver `ESRI Shapefile'
Simple feature collection with 10 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84

plot(suburbs$geometry)

tail(suburbs)

Simple feature collection with 6 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.46344 ymin: 55.63174 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
   id bydel_nr                      navn areal_m2                       geometry
5   6        6                   Vanløse  6699011 MULTIPOLYGON (((12.4982 55....
6   5        5                     Valby  9234177 MULTIPOLYGON (((12.52434 55...
7   4        4 Vesterbro-Kongens Enghave  8472602 MULTIPOLYGON (((12.54553 55...
8   1        1                  Indre By 10488466 MULTIPOLYGON (((12.72897 55...
9   9        9                Amager Øst  9820410 MULTIPOLYGON (((12.63082 55...
10  2        2                  Østerbro  9858740 MULTIPOLYGON (((12.59777 55...

# write_rds(suburbs, 'data/CPHsuburbs.rds')

suburbs$id

 [1]  3  7  8 10  6  5  4  1  9  2

Attribute data

Next let’s bring in the attribute data. I read the data from a google sheet where my colleagues and I can edit it. You can load it from there if you have a googlesheets4 package installed, or you can use the read_csv() function to read the wc.csv provided in the data folder

# Uncomment the lines below to read data from GDrive wc <- read_sheet('https://docs.google.com/spreadsheets/d/1iFvycp6M6bF8GBkGjA2Yde2yCIhiy5_slAkGF-RUF7w/edit#gid=0', col_types = 'cnnnnnnnn')
# write_csv(wc, 'data/wc.csv')
wc <- read_csv("data/wc.csv")
wc

# A tibble: 100 x 9
   Suburb          suburb_id flats wc_access  bath  year bath_communal_ct wc_communal_ct hot_water
   <chr>               <dbl> <dbl>     <dbl> <dbl> <dbl>            <dbl>          <dbl>     <dbl>
 1 Indre by                1 16280     11310  3800  1950               NA             NA      3390
 2 Christianshavn          1  5490      3900   900  1950               NA             NA       850
 3 Voldkvartererne         1 13460     12690  4560  1950               NA             NA      4290
 4 Østerbro                2 30820     28900 13750  1950               NA             NA     10750
 5 Indre Nørrebro          3 28700     24380  5910  1950               NA             NA      4010
 6 Ydre Nørrebro           3 21710     20410  5800  1950               NA             NA      4800
 7 Vesterbro               4 25850     23930  3730  1950               NA             NA      2480
 8 Kongens Enghave         4  6270      6240  4240  1950               NA             NA      3630
 9 Valby                   5 14430     13970  8190  1950               NA             NA      6980
10 Vigerslev               5  7700      7580  5050  1950               NA             NA      3850
# ... with 90 more rows

Spatial resolution adjustment - data aggregation

Data on access to hygienic facilities and other water resources in Copenhagen now looks good and tidy, but its spatial resolution is higher than the provided polygons (as in we have multiple rows that all fit within one suburb id). We therefore use the group_by() function to aggregate the data by id before we continue with any spatial operations. Given that the dataset is in fact a time-series, and each kvarter has a record for a given year or decade, we need to group first by the year and then only by id.

While aggregating the finer scale data into larger units, it is convenient to generate some statistics, such as percentages of flats that have bath and wc and hot water access within each suburb. We do this using the summarize() function below.

wcdata <- wc %>% group_by(year, suburb_id) %>% summarize(flats = sum(flats), bath = sum(bath), pct_bath = bath/flats * 100, wc_access = sum(wc_access), pct_wc = wc_access/flats * 100, warmH20 = sum(hot_water), 
    pct_wH20 = warmH20/flats * 100, communal_wc = sum(wc_communal_ct), communal_bath = sum(bath_communal_ct))
wcdata

# A tibble: 50 x 11
# Groups:   year [5]
    year suburb_id flats  bath pct_bath wc_access pct_wc warmH20 pct_wH20 communal_wc communal_bath
   <dbl>     <dbl> <dbl> <dbl>    <dbl>     <dbl>  <dbl>   <dbl>    <dbl>       <dbl>         <dbl>
 1  1950         1 35230  9260     26.3     27900   79.2    8530     24.2          NA            NA
 2  1950         2 30820 13750     44.6     28900   93.8   10750     34.9          NA            NA
 3  1950         3 50410 11710     23.2     44790   88.9    8810     17.5          NA            NA
 4  1950         4 32120  7970     24.8     30170   93.9    6110     19.0          NA            NA
 5  1950         5 22130 13240     59.8     21550   97.4   10830     48.9          NA            NA
 6  1950         6 10260  6780     66.1     10120   98.6    6270     61.1          NA            NA
 7  1950         7 27260 14790     54.3     26770   98.2   13280     48.7          NA            NA
 8  1950         8 19270 15000     77.8     18980   98.5   14690     76.2          NA            NA
 9  1950         9 23960 12470     52.0     22950   95.8   11210     46.8          NA            NA
10  1950        10 18000  9030     50.2     16200   90      7800     43.3          NA            NA
# ... with 40 more rows

# write_rds(wcdata, 'data/CPH_wcdata.rds')

Join the aggregated attribute data with its spatial representations

Now we can join the data with the spatial polygons

wc_spatial <- suburbs %>% merge(wcdata, by.x = "id", by.y = "suburb_id")
wc_spatial

Simple feature collection with 50 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
First 10 features:
   id bydel_nr     navn areal_m2 year flats  bath pct_bath wc_access   pct_wc warmH20 pct_wH20 communal_wc communal_bath                       geometry
1   1        1 Indre By 10488466 1981 26413 14035 53.13671     22546 85.35948      NA       NA          NA            NA MULTIPOLYGON (((12.72897 55...
2   1        1 Indre By 10488466 1950 35230  9260 26.28442     27900 79.19387    8530 24.21232          NA            NA MULTIPOLYGON (((12.72897 55...
3   1        1 Indre By 10488466 1965 32470 12780 39.35941     32450 99.93840      NA       NA        9510          2280 MULTIPOLYGON (((12.72897 55...
4   1        1 Indre By 10488466 1970 30440 11386 37.40473     22381 73.52497      NA       NA          NA            NA MULTIPOLYGON (((12.72897 55...
5   1        1 Indre By 10488466 1955 33490  9740 29.08331     33430 99.82084    9130 27.26187          NA            NA MULTIPOLYGON (((12.72897 55...
6   2        2 Østerbro  9858740 1950 30820 13750 44.61389     28900 93.77028   10750 34.87995          NA            NA MULTIPOLYGON (((12.59777 55...
7   2        2 Østerbro  9858740 1955 30540 14060 46.03798     30510 99.90177   11660 38.17944          NA            NA MULTIPOLYGON (((12.59777 55...
8   2        2 Østerbro  9858740 1965 32210 17340 53.83421     32190 99.93791      NA       NA        4670           920 MULTIPOLYGON (((12.59777 55...
9   2        2 Østerbro  9858740 1970 32464 17728 54.60818     28139 86.67755      NA       NA          NA            NA MULTIPOLYGON (((12.59777 55...
10  2        2 Østerbro  9858740 1981 32448 20255 62.42295     29343 90.43084      NA       NA          NA            NA MULTIPOLYGON (((12.59777 55...

Now that we have a merged spatial dataset with attributes, let’s review what attributes are available for visualisation

# Review the column names to see what new columns you have created
names(wc_spatial)

 [1] "id"            "bydel_nr"      "navn"          "areal_m2"      "year"          "flats"         "bath"          "pct_bath"      "wc_access"     "pct_wc"        "warmH20"       "pct_wH20"     
[13] "communal_wc"   "communal_bath" "geometry"

There is the suburb polygon data, such as id, bydel_nr, navn and areal_m2, and there is also the attribute data such as year, flats, bath ,etc. This gives us lots of choices for display. Lets put the data in a map.

Plot the data on the map

Let’s start by plotting one year alone, to learn how the map works.

Flats and water resources in 1950

Run the whole chunk below, and once it renders, look at the map. Afterwards, try changing the definition of what is to be displayed on line 116. For example, replace "flats" for some other column, such as "pct_bath", or "wc_access" to see how the map changes. To modify the legend, you can modify line 118 where we describe style. Replace style = "jenks" with "pretty", or "equal", or "quantile". What happens to your classification?

wc1950 <- wc_spatial %>% filter(year == 1950)

library(tmap)
tmap_mode(mode = "plot")
tm_shape(wc1950) + tm_borders(col = "black", lwd = 1) + tm_polygons("flats", id = "navn", style = "jenks") + tm_legend(legend.position = c("RIGHT", "TOP")) + tm_compass(position = c("RIGHT", "BOTTOM"), 
    type = "rose", size = 2) + tm_scale_bar(position = c("RIGHT", "BOTTOM"), breaks = c(0, 2, 4), text.size = 1) + tm_credits(position = c("RIGHT", "BOTTOM"), text = "Adela Sobotkova, 2021") + tm_layout(main.title = "Copenhagen 1950 situation", 
    legend.outside = FALSE)

Flats through time

Now, that you have mastered visualization of a single year, let’s plot all the years we have available!

#Check what years we have
unique(sort(wc_spatial$year))

[1] 1950 1955 1965 1970 1981

tmap_options(limits = c(facets.view = 5)) # we want to view 5 periods
tmap_mode(mode = "view" )
tm_shape(wc_spatial) +
  tm_facets(by = "year",
           ncol=5)+
  tm_polygons("flats",
              id = "navn",
             style = "jenks")+
  tm_layout(main.title = "Copenhagen Flats",
            legend.outside = TRUE)

Flats per square kilometer

Now that we have a spatial object, we can create new columns, for example utilizing the shape area to calculate the density of flats per sq km.

wc_spatial <- wc_spatial %>% mutate(area_km2 = areal_m2/1e+06, flat_per_km = flats/area_km2)

And now we can view these new columns in the map

tmap_mode(mode = "view")
tm_shape(wc_spatial) + tm_facets(by = "year", ncol = 5) + tm_polygons("flat_per_km", n = 5, style = "jenks")  #+

Access to private baths: percentages of households that have it per suburb

Lets calculate the baths and toilets available per square kilometer per each suburb

tmap_mode(mode = "view")
tm_shape(wc_spatial) + tm_facets(by = "year", ncol = 5) + tm_polygons("pct_bath", id = "navn", style = "pretty", title = "% of flats with <br> access to bath")

Access to private toilets: percentages of households that have it per suburb

library(tmap)
tmap_options(limits = c(facets.view = 5))  # we want to view 5 periods
tmap_mode(mode = "view")
tm_shape(wc_spatial) + tm_facets(by = "year", ncol = 5) + tm_polygons("pct_wc", id = "navn", style = "pretty", title = "% of flats with <br>access to WC")

Follow the example and develop your own analysis and map

You can also calculate the number of baths and toilets per square kilometer and plot in the map

wc_spatial <- wc_spatial %>% mutate(bath_per_km = bath/area_km2, wc_per_km = wc_access/area_km2)

..or continue with communal resources and warm water

Why not practice and try plotting the flats that have access to communal baths and wc, and or hot water? Create your own map here, following the examples above.